AD-A054  929 


UNCLASSIFIED 


AEROSPACE  CORP  EL  SE6UND0  CALIF  IVAN  A GETTING  LABS  F/G  20/5 

SPIKES  A COMPUTER  MODEL  FOR  THE  H2(D2)  ♦ F2  PULSED  CHEMICAL  LAS— ETC(U) 
APR  78  J J HOUGH  F04701-77-C-0078 


TR-0078< 3940-01) -4 


SAMSO-TR-78-84 


' AD  No. 

ODC  FILE  COPY  AD  A 0 5 49  29 


SPIKE:  A Computer  Model  for  the  H2(D2)  + F2 

Pulsed  Chemical  Laser 


J.  J.  T.  HOUGH 

Aerophysics  Laboratory 
The  Ivan  A.  Getting  Laboratories 
iThe  Aerospace  Corporation 
\ El  Segundo,  Calif.  90245 


14  April  1978 


Interim  Report 


APPROVED  FOR  PUBLIC  RELEASE; 
DISTRIBUTION  UNLIMITED 


Prepared  for 

SPACE  AND  MISSILE  SYSTEMS  ORGANIZATION 
AIR  FORCE  SYSTEMS  COMMAND 
Los  Angeles  Air  Force  Station 
P.O.  Box  92960,  Worldway  Postal  Center 


This  interim  report  was  submitted  by  The  Aerospace  Corporation, 

El  Segundo,  CA  90245,  under  Contract  No.  F04701-77-C-0078  with  the  Space 
and  Missile  Systems  Organization,  Deputy  for  Advanced  Space  Programs, 

P.O.  Box  92960,  Worldway  Postal  Center,  Los  Angeles,  CA  90009.  It  was 
reviewed  and  approved  for  The  Aerospace  Corporation  by  W.  R.  Warren,  Jr., 
Director,  Aerophysics  Laboratory.  Lieutenant  Dara  Batki,  SAMSO/YCPT, 
was  the  project  officer  for  Advanced  Space  Programs. 

This  report  has  been  reviewed  by  the  Information  Office  (OI)  and  is 
releasable  to  the  National  Technical  Information  Service  (NTIS).  At  NTIS, 
it  will  be  available  to  the  general  public,  including  foreign  nations. 

This  technical  report  has  been  reviewed  and  is  approved  for  publication. 
Publication  of  this  report  does  not  constitute  Air  Force  approval  of  the  report's 
findings  or  conclusions.  It  is  published  only  for  the  exchange  and  stimulation 
of  ideas. 


Dara  Batki, 
Project  Officer 


3 

2nd  Lt, 


USAF 


Robert  W.  Lindemuth,  Lt  Col,  USAF 


Chief,  Technology  Plans  Division 


r 


FOR  THE  COMMANDER 


FLOYD  R.  STUART.  Col,  USAF 
Deputy  for  Advanced  Space 
Programs 


TITLE  TO  33  VTRBALIMSD  lJ;  DJ/J-TA. 

FIELD  6 AND/OR  7 HOT  PUT  IKTO  SYSTEM  AT  S TAT TOM  T-3. 

Ol^a  G.  /DDE -TOD 

G2b  TITLE  VERBALIZED  BY  DDC-TA. 


(Mote  to  bo  removed  from  document  vhen  typist  at  T-l) 
puts  Fields  6 end/ or  7 into  system). 


UNCLASSIFIED  


SECURITY  CLARIFICATION  OF  THIS  PAGE  (•tt'h.n  Data  Entarad) 


EPORT  DOCUMENTATION  PAGE 


READ  INSTRUCTIONS 
BEFORE  COMPLETING  FORM 


2.  GOVT  ACCESSION  NO.I  3.  RECIPIENT'S  CATALOG  NUMBER 


[KE:  A COMPUTER  MODEL  FOR  THE 
(E^)  + F^  PULSED  £ HE  MIC  A L J^ASER. 


Interim 


1 


7.  author^; 


V Joseph  J. 


9.  PERFORMING  ORGANIZATION  NAME  ANO  ADDRESS 

The  Aerospace  Corporation 
El  Segundo,  Calif.  90245 


10.  PROGRAM  ELEMENT.  PROJECT,  TASK 
ARg-IL*  WBRII  UIIIP  MUMPERS 


78  ! ✓ 


It.  CONTROLLING  OFFICE  NAME  AND  ADDRESS  f 

Space  and  Missile  Systems  Organization  / / J 
Air  Force  Systems  Command 

Los  Angeles,  Calif.  90009  


14.  MONITORING  AGENCY  NAME  & ADDRESS^//  different  from  Controlling  Office)  15.  SECURITY  CLASS,  (of  this  report) 

Unclassified 

1 5a.  DECLASSIFICATION/ DOWN  GRADING 
SCHEDULE 


16.  DISTRIBUTION  STATEMENT  (of  thla  Raport) 


Approved  for  public  release;  distribution  unlimited 


17.  DISTRIBUTION  STATEMENT  (of  the  abstract  entered  in  Block  20,  if  different  from  Report) 


19.  KEY  WORDS  (Continue  on  reverse  aide  If  necessary  and  identify  by  block  number) 

Lasers  Pulsed  Lasers 

Chemical  Lasers  Infrared  Lasers 

HF 

Computer  Model 


20.^  ABSTRACT  (Continue  on  reverse  aide  if  necessary  and  identify  by  block  number) 

A^comprehensive  rate  equation  model  for  a pulsed  laser  pumped  by  the 
H2(E>2)-F2  reaction  is  described.  This  computer  simulation  incorporates 
up  to  216  kinetic  reactions  to  describe  the  laser  medium,  which  is  assumed 
to  be  homogeneous  and  contained  in  a Fabry-Perot  optical  cavity.  The  first 
fifteen  vibrational-rotational  P-branch  transitions  from  each  of  the  lowest 
ten  vibrational  bands  are  monitored  for  possible  laser  action.  The  distinctive 
features  incorporated  in  this  model  include:  (1)  the  elimination  of  the  gain- 
equals-loss  assumption,  so  that  it  is  possible  to  predict  cavity  transients 


DD  F0RM  1473 

IPA  CSimilE  > 


UNCLASSIFIED 

SECURITY  CLASSIFICATION  OF  T>l I S PAGE  (Whan  Data  I 


:URI1Y  CLASSIFICATION  OF  TMI 

#0/ 977 


UNCLASSIFIED 

SECURITY  CLASSIFICATION  OF  THIS  PAOEftWmi  Data  Bnlotod) 


19.  KEY  WORDS  (Continu'd) 


20.  ABSTRACT  ( Continued ) 

^and  to  observe  gain -intensity  interactions;  (2)  a detailed  model  of  the 
Lorentz  broadening,  including  V,  J dependence,  as  well  as  distinction 
between  various  perturbing  species;  (3)  multiquantum  HF(DF)  self- 
deactivation, which  in  the  HF  case  was  found  important  to  achieve  good 
agreement  between  theory  and  experiment. 

Comparisons  of  the  prediction  of  this  model  with  several  small-scale  experi- 
ments have  indicated  generally  good  agreement  (typically  within  ~20%). 

A 

\ 


UNCLASSIFIED 


SECURITY  CLASSIFICATION  OF  THIS  RAOtfWAan  Dmtm  Enfrod) 


PREFACE 
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at  Michigan  State  University  with  R.  L.  Kerber  and  was  supported  by  the 
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with  experiments,  was  carried  out  at  The  Aerospace  Corporation  and  was 
supported  by  the  Defense  Advanced  Research  Projects  Agency,  the  Naval 
Sea  Systems  Command,  and  the  Air  Force  WeaponsLaboratory  under 
U.S.  Air  Force  Space  and  Missile  Systems  Organization  (SAMSO)  Contract 
No.  F04701 -77-C -0078. 

This  code  originally  was  not  intended  for  distribution;  hence,  it  may 
not  contain  some  convenience  measures  that  might  be  desirable  from  the 
standpoint  of  the  unfamiliar  user.  Sufficient  detail  is  provided  in  this  report, 
however,  such  that,  with  some  effort,  the  user  will  be  able  to  modify  the 
program  for  his  specific  purposes. 
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I.  INTRODUCTION 

/ \ 

/ SPIKE,  a comprehensive  computer  model  of  the  HF(DF)  pulsed  chemical 
laser  pump^d'by  the  H2(D2)  + F^  chain  reaction,  is  described.  This  model 
incorpcrfates  up  to  216  kinetic  reactions  and  takes  into  detailed  account  the 
effects  of  pressure  line -broadening  based  on  the  results  of  a recent  survey 
of  the  literature  (Ref.  1).  The  method  of  computation  utilized  in  this  model 
did  not  require  the  gain-equals -loss  assumption  and  is  capable  of  predicting 
transients  and  relaxation  oscillations  within  the  laser  cavity  (Ref.  2). 

Recent  comparisons  of  the  predictions  of  this  model  with  several 
small-scale  pulsed  H2~F2  laser  (~25-cm  gain  length)  experiments  revealed 
generally  good  agreement  (Ref.  2).  Comparisons  with  larger  scale  experi- 
ments (~100-cm  gain  length),  however,  have  not  shown  as  good  agreement 
in  the  prediction  of  absolute  energy,  although  relative  performance  with 
parameter  variations  still  showed  comparable  results  (Ref.  3).  This  lack 
of  agreement  is  indicative  of  possible  deficiencies  in  the  model.  It  has  been 
suggested  that  such  phenomena  as  photoionization,  photodissociation,  and 
medium  nonuniformity,  which  have  been  neglected  in  typical  modeling 
calculations,  may  have  significant  impact  on  the  larger  scale  lasers.  These 
possibilities  are  still  under  investigation. 

Several  recent  requests  for  SPIKE,  reflect  the  need  for  this  "users 
guide."  The  objectives  in  this  report  are: 

1.  To  explain  the  theories  from  which  the  calculations 
executed  in  the  code  are  formulated. 

2.  To  present  in  systematic  form  the  formal  and  logical 
structures  of  the  code. 

3.  To  describe  the  procedures  necessary  in  order  to  obtain 
information  from  this  code  with  minimal  lag  time. 

4.  To  provide  sufficient  detail  such  that,  with  some  prepara- 
tion, modifications  to  the  code  can  be  made  so  that  it  can 
be  tailored  to  particular  needs. 
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II.  MODEL  FORMULATION 

The  formulation  of  the  chemical  laser  computer  model  is  described  in 
this  section.  The  reactions  used  to  represent  the  chemical  processes  are: 

1 . H2  + F^  chain 

F + H2^HF(v)  + H 
H + F2^HF(v)  + F 

2.  Vibrational  -translational  (VT ) deactivation 

HF  (v)  + M^HF(v')  + M 
H2(v)  + M=^H2(v  - 1)  + M 

3.  Vibrational -vibrational  (VV)  quantum  exchange 

HF(v)  +HF(v')^HF(v  + 1)  +HF(v'  - 1) 

HF(v)  +H2(v')^HF(v  + 1)  +H2(v'  - 1) 

4.  Dissociation-recombination 

F 2 + M ^M  + F + F 

H2+M^M+H+H 
HF  + M^M  + H + F 
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The  major  features  in  the  models  are: 

1.  The  dominant  kinetic  processes  are  represented  by  the 
reaction  system  suggested  by  Cohen  (Ref.  4)  (Table  A-l). 

2.  The  reacting  mixture  is  homogeneous  and  is  contained  in 
a Fabry-Perot  laser  cavity. 

3.  The  rotational  populations  for  each  vibrational  level  have 
a Boltzmann  distribution  at  the  translational  temperature. 

4.  All  possible  transitions  within  a band  are  assumed  to  have 
low  initial  intensities  that  grow  if  the  gain  rises  above 
threshold.  Lasing  is  always  assumed  to  be  in  the  P-branch. 
Initial  intensity  levels  can  be  selected  individually  or  set 
proportional  to  the  spontaneous  emission  rate. 

5.  Initiation  is  modeled  by  the  introauction  of  a finite  concen- 
tration of  F atoms  into  the  gas  mixture,  or  through  the  use 
of  a flash-lamp  option. 

The  chemical  reactions  are  written 


Eo  .N. 
ri  x 


-r  i 


(1) 


where  N.  is  the  molar  concentration  of  species  i,  a . and  3 . are  stoichio- 
x ri  ri 

metric  coefficients,  and  k and  k are  forward  and  backward  rate  coeffi- 

r -r 

cients.  The  rate  of  change  of  concentration  for  nonlasing  molecules  is 


(2a) 


and,  for  HF  molecules, 


dN. 


= Xi  + Xrad<v'  J>  - W ' *•  JL 


) + A(v,  J) 


(2b) 


1* 


where  the  terms  are  rates  of  change  in  concentration  as  a result  of 

lasing  into  and  out  of  level  (v,  J).  The  lower-level  rotational  quantum  numbers 

are  J and  JT  for  the  transitions  v + 1 — > v and  v — * v - 1,  respectively.  The 
-Ll 

net  rate  of  spontaneous  emission  into  level  (v,  J)  is  given  by  A(v,  J).  The 
chemical  reactions  yield  a concentration  change 


(3) 


and 

Xrad(v’  J)  = g(v’  J)  f(v’  J)  (4) 

where  g(v,  J)  is  the  gain  on  the  v + 1 — > v transition  with  lower  level  J and 
f(v,  J)  is  the  lasing  flux  on  the  same  transition.  The  rate  equation  for  the 
lasing  flux  is 

^ dt  J)  = C(J)  [g(v’  J)  - athrlf(v’  J)  (5) 

where  c is  the  speed  of  light  and 


O',  = - ■=■=■  ln(R  RT  ) 

thr  2L,  o L 


(6) 


where  L is  the  length  of  the  active  medium,  i is  the  mirror  spacing,  and 
Rq  and  R^  are  the  mirror  reflectivities.  The  gain  (Ref.  5)  is 


hN 


g(v,  J)  = i»c(v,  J)  J)  B(v,  J) 

( ZJ  + l\j^  _ 

* y2 j - ly iNHF(v+l,  J-l) 


HF (v,  J) 


(7) 
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where  u)c(v,  J)  is  the  wave  number  of  the  transition,  B(v,  J)  is  the  Einstein 
isotropic  absorption  coefficient  based  on  the  intensity,  and  0(v,  J)  is  the 
Voigt  profile  at  line  center  (Ref.  6).  The  first  term  of  Eq.  (5)  determines 
the  rate  of  increase  in  the  intensity  of  the  radiation  field  within  the  laser 
cavity;  the  second  term  gives  the  rate  energy  that  is  lost  from  the  cavity.  The 
lost  energy  includes  that  extracted  through  the  output  coupler  as  laser  output, 
as  well  as  real  losses  resulting  from  such  mechanisms  as  absorption, 
scattering,  and  extraneous  reflections. 


A Boltzmann  distribution  at  the  translational  temperature  (T)  is  taken 
for  the  rotational  popul  ations ; hence, 


_ _ _ _ /2 J + 1 

HF(v,  J)  = NHF(v)|Qv(T)^|e 


-hcEj/kT 


(8) 


where  the  values  of  the  rotational  partition  function  Q^(T)  and  the  rotational 
energy  Ej  are  computed  from  Mann  et  al.(Ref.  7). 

The  energy  equation  for  a constant  density  gas  is 


dT 

dt 


N.C 
1 v. 

l 


-P, 


dN. 
dt  i 


(9) 


where  C is  the  molar  specific  heat  at  constant  volume,  H.  is  the  molar 
v.  1 

enthalpy  of  species  i,  and  P is  the  output  lasing  power  per  unit  volume. 


The  output  power  in  the  v + 1 — » v band  is 


where  the  only  cavity  loss  is  assumed  to  be  the  laser  output.  In  making 
comparisons  with  experimental  measurements,  however,  real  losses  must 
be  accounted  for.  One  way  to  do  this  is  to  assume  that  fraction  of  the  output 
attributable  to  the  "total  reflector"  RT  (as  distinguished  from  the  output 
coupler  Rq)  corresponds  to  cavity  losses.  The  transmissivity  of  the  mirror 
(R  ) may  be  adjusted  to  approximate  cavity  losses  resulting  from  scattering, 
transmission,  and  absorption  at  various  optical  surfaces.  The  useful  frac- 
tion of  the  total  output  power  is  then  (Ref.  6) 


Ptotal^-V 

useful  [1  + (Rq/R1j)1 /2]  [ 1 - (Ro/RL)1/2j 


From  the  numerical  integration  of  Eqs.  (2),  (5),  and  (9)  by  the  modified 
Adams -Moulten  method  of  Gear  (Ref.  8),  the  time  evolution  of  the  species 
concentrations,  temperature,  pressure,  the  gain  on  all  transitions,  and  the 
intensities  on  all  lasing  transitions  are  determined.  The  laser  energy  ex- 
tracted in  each  band  is  then  determined  by  integrating  the  power 


rt 

=/c 


Pr  dt 
Lv 


(12a) 


where  t is  the  length  of  laser  pulse  and  the  total  pulse  energy  is 


=I>V 


(12b) 


J 


T 
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III.  CODE  DESCRIPTION 


The  broad  features  of  SPIKE  are  described  in  this  section.  The  code 
is  written  in  FORTRAN.  A complete  listing  of  the  main  program  and  the 
subroutines  is  given  in  Appendix  C.  In  addition  to  the  explanations  provided 
in  this  and  the  subsequent  sections,  extensive  comment  statements  are  in- 
cluded throughout  the  body  of  the  listing  to  provide  further  clarification. 

The  program  consists  of  a main  (or  controlling),  program,  i.e.,  SPIKE, 
and  six  subroutines:  INPUT  1,  LORENTZ,  GAIN.  DIFSUB,  DIFFUN,  and 
RATE.  SPIKE  obtains  input  data  by  calling  subroutine  INPUT  1 and  initiates 
such  computational  variables  as  concentration,  photon  flux,  temperature 
and  time;  defines  the  cavity  conditions;  obtains  reaction  rate  constants  and 
necessary  spectroscopic  and  thermodynamic  data  from  data  files;  updates 
the  pressure  line -broadening  parameters  through  the  use  of  the  subroutine 
LORENTZ;  and  provides  the  program  control  parameters  that  direct  the 
program  flow,  as  well  as  the  integration  process.  Subroutine  (function) 

GAIN  is  used  to  compute  the  medium  gain  of  the  laser  for  any  given  allowed 
vibrational -rotational  P-branch  transition  by  means  of  Eq.  (7).  The  photon 
fluxes  corresponding  to  the  seven  transitions  within  each  band  with  the 
highest  gain  are  monitored  for  possible  laser  action. 

Time  progression  within  the  model  proceeds  along  the  integration 
steps.  Numerical  integration  is  accomplished  through  the  use  of  the  sub- 
routine DIFSUB,  which  utilizes  the  modified  Adams  - Moulten  technique 
presented  by  Gear  (Ref.  8).  Derivatives  required  for  the  integration  are 
computed  from  Eqs.  (2),  (5),  and  (8)  and  obtained  by  means  of  subroutine 
DIFFUN.  Laser  initiation  by  electrical  discharge  (or  flash  photolysis)  may 
be  modeled  through  a flash-lamp  option  contained  in  the  main  program, 
which  provides  a time -dependent  F-atom  production  rate  based  on  the  input 
power  (intensity)  profile. 
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Substantial  savings  in  computation  time  and  computer  core  storage 
requirements  were  obtained  with  the  present  code  formulation,  wherein  the 
rate  equations  were  input  as  an  integral  part  of  the  program,  and  the  rate 
coefficients  were  precomputed  and  stored  on  disk  for  retrieval  during  pro- 
gram execution.  In  the  event  that  minor  modifications  to  the  rates  are 
desired,  rather  than  re-creating  the  entire  rate  file,  these  minor  changes 
are  accomplished  through  the  use  of  subroutine  RATE. 

Figure  1 is  a schematic  flow  diagram  of  the  computer  program.  This 
figure  is  an  overview  of  the  chain  of  commands  and  decisions  used  to  obtain 
the  integrated  values  of  the  species  concentrations  within  the  laser  medium, 
as  well  as  the  laser  power  density  and  temperature  as  a function  of  time. 

The  variables  and  symbols  used  in  Fig.  1 and  in  the  computer  program 
are  defined  in  Tables  1 and  2.  Bookkeeping  limitations  and  computational 
practicability  demand  that  the  many  time -dependent  variables  be  grouped  in 
an  array.  In  addition,  since  the  dimension  of  an  array  is  changeable,  the 
model  has  the  flexibility  to  alter  the  number  of  variables  under  consideration. 

The  array  is  two  dimensional  and  is  designated  Y(I,  J),  J = 1 N,  where 

Y(l,  J)  represents  the  variable  in  question,  N is  the  total  number  of  such  vari- 
ables, and  Y(I,  J)  and  I > 1,  are  related  to  the  (I  - l)th  derivatives  of  Y(l,  J). 
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Table  1.  Identification  of  Variables  Y(l,  N)  in 
H2(D2)  + F2  Model 


Y(  1 , N) 


1-9 

10-11 

12-14 

15-16 

17 

18-20 

21 

22 

23 

24-79 


_ 3 

HF(v)[DF(v)l  concentration,  mol-cm  , v = N-l 

_3 

DF(v)  concentration,  mol-cm  , v = N-l 

_3 

N2(v)  concentration,  mol-cm  , v = N-12 

Undefined,  reserved  for  additional  variables 

-3 


H(D)  - atom  concentration,  mol-cm 


v = N-  18 


_3 

H2(v)  [D2(v)]  concentration,  mol-cm 

-3 

F -atom  concentration,  mol-cm 

_3 

F2  concentration,  mol-cm 
Translational  temperature,  K 

Photon  flux,  f(v,  J),  mol-cm  -sec  , v=  | (N - 17 )/7  | , 

J = Jv  - 4,  ....  Jv  +2,  where  JV  indicates  the 
max  max  max 


transition  of  maximum  gain. 
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Table  2.  Nomenclature 


Symbol  in 
Text 

Symbol  in 

Computer  Program 

Definition 

Afv,  J) 

A( V.  J) 

Einstein  isotropic  coefficient  for 
spontaneous  emission.  1 /molecule- 
sec 

B( v,  J) 

B( V.  J) 

Einstein  isotropic  intensity  absorption 
coefficient,  cm^/molecule  - J - sec 

c 

C 

Speed  of  light.  2.  997925  x 10*^  cm/sec 

cv. 

1 

cvi 

Molar  specific  heat  at  constant  volume 
of  species  i.  cal/mol-K 

*5 

E(V.  J) 

Rotational  energy  of  state  V.  J.  cm  * 

f(v.  .n 

FLUX(V.  .F) 

Photon  flux,  mol -cm  ^-sec  * 

K(v.  I) 

ALPH A( V.  .1) 

Gain  of  transition  (v  ♦ 1.  .1  - l ) 

— * (v.  J ).  cm  ‘ 1 

h 

- 

Planck's  constant,  6.625b  x 10"^  J-sec 

h 

Specific  enthalpy,  kcal/g 

H. 

1 

EHPY! 

Molar  enthalpy  of  species  i.  kcal/mol 

I 

CRRNT 

Discharge  current.  A 

k 

K 

-23 

Boltzmann  s constant.  1 3H054  x 10 

F-K-1 

k . k 
r - r 

l\l- 'll.  KltU 

f orward  and  ba<  kward  rati-  constants, 
in  terms  of  moles,  centimeters,  and 
seconds 

L 

LNTH 

Length  of  active  medium,  cm 

N. 

l 

Y(  1 . I) 

Concentration  of  species  i.  mol-cm  ^ 

dN.  /dt 

i 

DERV1YH) 

Time  derivative  of  N..  mol-cm  ^ sec  * 

na 

NA 

Avogadro's  number.  6.  02252  x 10^ 

molecules -mol " ^ 

PL 

TOWER 

Power  density  of  laser  output.  W-ctn  ^ 

P. 

PIN 

Power  input  from  initiation.  W-cm  ^ 

0 

1 < 

Q(V.  T) 

Rotational  partition  function  for  level  v 

rorl 

RO.  RL 

Mirror  reflectivities 

R 

R 

Universal  eas  constant.  1 . 98725  cal  - 
mol-l-K'f 

t 

T 

Time,  sec 

T 

Y(  1.  23) 

Temperature.  K 

a ..  P . 

ri  ri 

- 

Stoichiometric  coefficients  of  reaction  r 

0,thr 

THGAIN 

Threshold  gain,  cm  * 

0(V.  J) 

PHI(V.  J) 

Normalized  line  profile  of  transition 
(v  f 1.  J - 1 ) — * (v.  J ).  cm 

tr 

TAU( V.  J) 

Rotational  relaxation  time  constant,  sec 

uu  (V.  J) 
c 

WC(V.  J) 

Wave  number  of  transition  (v  ♦ 1.  .1  IF 

— * |v.  J).  cm”  * 

I 


IV.  BASIC  INSTRUCTIONS  FOR  USING  SPIKE 


The  common  input  preparations  needed  for  most  runs  are  discussed  in 
this  section.  The  input  parameters  are  classified  as  follows:  (1)  initial  gas 
conditions,  (2)  optical  cavity  parameters,  (3)  output  control,  (4)  integration 
control  parameters,  and  (5)  initiation  mechanism. 

A sample  set  of  input  data  cards  is  shown  in  Fig.  2.  The  proper  loca- 
tion for  these  cards  is  immediately  after  the  end-of -record  card  that  follows 
the  program  dec  s.  This  format  corresponds  to  that  of  NAMELIST.  A 
summary  of  the  input  variables  is  given  in  Table  3. 

The  input  data  are  divided  into  seven  groups;  each  group  has  an 
identifying  name.  The  data  stream  for  each  group  is  initiated  by  the  sign  $ 
in  Column  2,  followed  by  the  group  name  (no  blank  must  appear  between  $ 
and  the  group  name)  and  then  the  data.  The  data  must  be  separated  by 
commas,  and  terminated  by  another  $ . However,  the  data  within  each  group 
can  be  listed  in  any  order  and  may  extend  to  more  than  one  record.  It 
should  be  emphasized  that  no  data  can  appear  in  Column  1 of  any  one  record, 
as  it  is  ignored  by  the  NAMELIST  READ  operation.  The  seven  data  streams 
are  described  below. 

A.  INITIAL  GAS  CONDITIONS 

The  initial  concentrations  of  the  species  in  the  gas  mixture  are  provided 
through  the  data  stream  SPECIES.  In  addition  to  the  reacting  species  F£, 

H^,  F,  and  H,  several  diluents  are  modeled,  i.e.,  He,  Ar,  N^,  SF^,  and  O^. 
These  concentrations  may  be  provided  in  any  units  (but  must  be  consistent) 
or  simply  as  relative  ratios.  An  initial  HF(0)  population  may  also  be  entered. 
Pressure  of  the  gas  mixture  in  Torr  (PRESS)  and  initial  gas  temperature  in 
degrees  kelvin  (TEMP)  are  input  by  the  data  stream  GAS. 


array . 


B. 


OPTICAL-CAVITY  PARAMETERS 


The  input  parameters  for  the  laser  optical  cavity  are  RO,  RL,  LNTH, 
LRATIO,  and  FLUX.  These  quantities,  contained  in  data  stream  CAVITY, 
are  defined  in  this  section. 

The  optical  cavity  is  assumed  to  be  Fabry-Perot,  with  mirror  reflec- 
tivities denoted  by  RO  and  RL  and  given  in  decimal  fractions.  RO  is  the 
output  coupler,  and  only  that  fraction  of  the  total  radiant  energy  that  passes 
through  this  mirror  is  included  in  the  calculation  of  laser  output.  Note  that 
the  "0"  in  the  symbol  RO  is  a "zero,  " not  the  letter  "O.  " The  gain  length  of 
the  active  medium  is  given  by  the  variable  LNTH  in  centimeters. 

LRATIO  is  the  ratio  of  the  gain  length  to  the  separation  of  the  mirrors, 
corresponding  to  the  L/f  factor  in  Eq.  (5).  This  parameter  accounts  for 
the  intervals  during  which  the  photons  traverse  regions  between  the  mirrors, 
which  are  neither  amplifying  nor  absorbing.  The  larger  the  "null"  region, 
the  fewer  passes  each  photon  makes  through  the  active  medium  during  the 
photon  residence  time  within  the  cavity.  This  phenomenon  reduces  the 
amplification  potential  of  the  cavity  and  can  result  in  appreciable  reduction 
in  laser  output.  It  has  been  observed  experimentally  at  Sandia  Laboratories.  * 

All  possible  transitions  within  a band  are  assumed  to  have  initial 

intensity  levels  that  grow  if  the  gain  rises  above  threshold.  'T'hese  initial 

intensity  levels  can  be  selected  individually  or  set  proportional  to  the 

spontaneous  emission  rate.  Our  standard  practice  has  been  to  set  all  initial 

- 14  2 

intensity  levels  at  10  mol/cm  -sec,  which  is  done  by  setting  the  variable 

FLUX  to  this  value.  Test  cases  where  this  value  was  varied  by  as  much 
3 

as  X 10  have  indicated  the  calculations  to  be  relatively  insensitive  to  this 
number  (Ref.  2).  The  value  of  FLUX  is  set  to  zero  to  make  small  signal 
gain  calculations.  (Alternatively,  the  threshold  gain  value  can  be  increased 
to  any  large  number  by  adjusting  RO,  RL,  and  LNTH.  ) 

*J.  B.  Moreno,  Sandia  Laboratories,  Albuquerque,  N.  M. , private 
communication,  April  1977. 
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c. 


OUTPUT  CONTROL 


TSKIP,  KBS,  and  IPEEK  are  three  input  parameters  that  control  the 
quantity  and  form  of  the  computer  output.  As  the  calculation  progresses, 
pertinent  data  are  printed  at  time  intervals  specified  by  the  parameter  TSKIP 
(seconds)  beginning  at  time  t = 0.  A final  print  is  also  obtained  at  pulse 
termination.  Contents  of  the  printed  output  are  controlled  by  the  parameter 
KBS,  which  may  have  values  of  1 < KBS  S 5.  The  larger  the  value  selected, 
the  more  print  is  obtained.  The  outputs  corresponding  to  the  various  values  of 
KBS  are  shown  in  Table  4. 


Table  4.  Definition  of  Parameter  KBS 


KBS 

Printed  Output 

1 

Time,  laser  power,  laser  energy,  gas  temperature,  species 
concentrations,  total  number  of  integration  steps  to  reach 
this  time,  and  the  last  step  size. 

2 

Output  of  KBS  = 1 plus  the  photon  flux,  power,  and  energy 
of  each  laser  transition. 

3 

Output  of  KBS  = 2 plus  the  gain  of  each  transition  at  this 
time . 

4 

Output  of  KBS  = 3 plus  the  self -broadening  coefficients 

Y(v,  J)  and  other  linewidth  data. 

5 

Output  of  KBS  = 4 plus  rate  coefficients  and  reaction  rates 
for  all  reactions. 

In  addition  to  printed  output,  an  array  (XXX)  is  written,  which  contains 
the  laser  power,  laser  energy,  time,  gas  temperature,  species  concentrations, 
and  gain  and  power  of  each  transition.  At  each  print  interval,  XXX  is  copied 
onto  Tape  13.  The  parameter  KKK  is  incremented  each  time  to  keep  track 
of  the  number  of  times  XXX  is  copied.  At  the  end  of  the  calculation.  Tape  13 
may  be  saved,  either  for  later  use  in  a plotting  program  or  for  a brief  scan 
of  the  calculated  results  (e.g.,  this  file  may  be  quickly  retrieved  and  studied 

If-  (’ 

> 


on  an  interactive  terminal  while  the  main  output  is  waiting  to  print).  If 
only  the  abbreviated  output  data  contained  in  XXX  is  needed,  the  main  print 
may  be  skipped  by  setting  IPEEK  = 1 (otherwise,  it  should  be  set  to  zero). 

The  resulting  output  will  be  two  prints  corresponding  to  the  initial  (t  = 0) 
and  final  (pulse  termination)  conditions. 

The  input  data  TSKIP.  KBS,  and  IPEEK  are  punched  in  data  stream 
OUTPUT1.  In  addition  to  these,  this  data  stream  also  contains  the  data 
TIMEL,  POFF,  and  TCHECK,  which  determine  the  point  at  which  computation 
is  stopped. 

Termination  of  the  calculation  occurs  with  the  onset  of  one  of  two 

events:  (1)  integration  reaches  an  imposed  time  limit  specified  by  TIMEL 

(seconds)  or  (2)  the  power  density  drops  below  some  preset  level,  given  by 
3 

POFF(W/cm  ).  In  the  latter  case,  to  avoid  terminating  the  calculation 
during  the  developing  phase  of  the  laser  power,  this  terminating  mechanism 
does  not  become  operative  until  an  interval,  TCHECK  (seconds),  has  elapsed. 
The  user  must  provide  values  for  TIMEL,  TCHECK,  and  POFF  as  part  of 
the  input  data.  In  the  case  of  a small -signal-gain  calculation,  where  power 
remains  at  zero,  termination  of  computations  is  controlled  only  by  TIMEL. 
The  value  of  TCHECK  should  be  set  larger  than  that  of  TIMEL. 

D.  INTEGRATION  CONTROL  PARAMETERS 

Numerical  integration  of  the  rate  equations  within  the  present  model  is 
accomplished  utilizing  a modified  Adams -Moulten  technique  presented  in 
Ref.  8.  The  corresponding  integration  subroutine  is  called  DIFSUB.  The 
reader  is  referred  to  the  original  reference  for  details  regarding  this  inte- 
gration package.  Under  ordinary  circumstances,  only  the  variables 
NVAR(  = N),  DELTA(^H),  HMIN,  H MAX,  EPS,  YMAX,  KFLAG,  and  JSTART 
will  be  of  concern.  NVAR  represents  the  number  of  integration  variables 
contained  in  the  program,  DELTA  is  the  step  size  to  be  attempted  on  the 
next  integration  step,  HMIN  and  HMAX  are,  respectively,  the  minimum  and 
maximum  step  sizes  that  will  be  used  for  the  integration.  EPS  is  an  error 
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test  constant,  which  represents  the  maximum  error  (estimate)  that  can  be 
tolerated  on  a step;  it  is  expressed  as  a "fraction"  of  the  maximum  value 
(YMAX)  of  each  integration  variable  (Y)  seen  up  to  this  point  in  the  integration. 
The  user  also  has  the  option  of  adjusting  the  values  of  YMAX  to  influence 
the  step  sizes.  KFLAG  is  an  output  flag  from  DIFSUB  that  gives  the  result 
of  the  previous  integration  attempt.  A KFLAG  value  of  +1  indicates  a success- 
ful step,  whereas  a negative  value  corresponds  to  an  unsuccessful  one. 

JSTART  defines  the  action  on  the  next  integration  attempt:  to  repeat  the 
previous  step,  to  continue  with  a new  step,  or  to  indicate  the  first  step.  A 
more  detailed  description  of  these  and  other  variables  is  given  in  Table  5; 
the  contents  of  which  are  taken  from  Ref.  8.  The  value  of  DELTA,  HMIN, 
HMAX,  EPS,  JSTART,  and  NVAR  must  be  provided  through  the  data  stream 


STEPS  to  make  a calculation. 


INITIATION 


Initiation  of  the  chemical  reactions  is  either  through  an  instantaneous 
introduction  of  F atoms  into  the  gas  mixture  or  a time -dependent  dissociation 
of  F^  molecules  by  a flash-lamp  option.  If  the  flash-lamp  option  is  selected, 
input  values  for  the  parameters  Zl,  Z2,  and  WPHOTO  must  be  provided, 
along  with  a flash-lamp  output  intensity  profile  (I  vs  t).  These  data  are 
entered  by  means  of  data  streams  FLASH1  and  FLASH2.  For  instantaneous 
initiation,  on  the  other  hand,  the  user  need  only  provide  an  initial  F -atom 
concentration  in  the  data  stream  SPECIES,  neglecting  the  data  streams 
associated  with  the  flash-lamp  option  FLASH1  and  FLASH2. 


The  formulation  of  the  flash -lamp  option  in  the  present  model  is  similar 
to  that  in  RESALE  (Ref.  6).  The  rate  of  dissociation  of  the  absorbing  species 


n&  is  given  by 


/dn  \ 

(ifja  = -exp(-Z2na)l 


A 


Table  5.  Definition  of  Parameters  in  DIFSUR 


I 


Parameter  Definition 

N = NVAR  The  number  of  first-order  differential  equations.  N may  be  decreased  on  later  calls  if  the 

number  of  active  equations  reduces,  but  it  must  not  be  increased  without  calling  with 
JSTART  = 0. 

T The  independent  variable. 

Y An  8 by  N array  containing  the  dependent  variables  and  their  scaled  derivatives.  Y(J  + 1,  I) 

contains  the  J-th  derivative  of  Y(I)  scaled  by  H /facto rial(J)  where  H is  the  current  step 
size.  Only  Y(l,  I)  need  be  provided  by  the  calling  program  on  the  first  entry. 

If  it  is  desired  to  interpolate  to  non  -mesh  points,  these  values  can  be  used.  If  the  current 
step  size  is  H and  the  value  at  T + E is  needed,  form  S = E/H,  and  then  compute 

NQ 

Y(I)  (T  f E)  = SUM  Y ( J + 1,  I)*S**J 
J = 0 

SAVE  A block  of  at  least  12*N  floating  point  locations  used  by  the  subroutines. 

H = DELTA  The  step  size  to  be  attempted  on  the  next  step.  H may  be  adjusted  up  or  down  by  the  program 

in  order  to  achieve  an  economical  integration.  However,  if  the  H provided  by  the  user  does 
not  cause  a larger  error  than  requested,  it  will  be  used.  To  save  computer  time,  the  user 
is  advised  to  use  a fairly  small  step  for  the  first  call.  It  will  be  automatically  increased 
later. 

HMIN  The  minimum  step  size  that  will  be  used  for  the  integration.  Note  that  on  starting  this  must 

be  much  smaller  than  the  average  H expected  since  a first-order  method  is  used  initially. 

UMAX  The  maximum  size  to  which  the  step  will  be  increased. 

EPS  The  error  test  constant.  Single-step  error  estimates  divided  by  YMAX(I)  must  be  less  than 

this  in  the  Euclidean  norm.  The  step  and/or  order  is  adjusted  to  achieve  this. 

MF  The  method  indicator.  The  following  are  allowed: 

0 An  Adams  predictor  corrector  is  used. 

1 A multi-step  method  suitable  for  stiff  systems  is  used.  It  will  also  work 
for  non-stiff  systems.  However,  the  user  must  provide  a subroutine 
PEDERV  which  evaluates  the  partial  derivatives  of  the  differential  equations 
with  respect  to  the  Y's.  This  is  done  by  call  PEDERV(T,  Y,  PW,  M).  PW  is 
an  N by  N array  which  must  be  set  to  the  partial  of  the  I-th  equation  with 
respect  to  the  J dependent  variable  in  PW(I,  J).  PW  is  actually  stored  in  an 

M by  M array  where  M is  the  value  of  N used  on  the  first  call  to  this  program. 

2 The  same  as  Case  1,  except  that  this  subroutine  computes  the  partial 
derivatives  by  numerical  differencing  of  the  derivatives . Hence  PEDERV  is 
not  called. 

YMAX  An  array  of  N locations  which  contains  the  maximum  of  each  Y seen  so  far.  It 

should  normally  be  set  to  1 in  each  component  before  the  first  entry  (See  the  de- 
scription of  EPS.  ) 

ERROR  An  array  of  N elements  which  contains  the  estimated  one-step  error  in  each 

component . 

KFLAG  A completion  code  with  the  following  meanings: 

*■  1 The  step  was  successful. 

-1  The  step  was  taken  with  H - HMIN.  but  the  requested  error  was  not  achieved. 

-2  The  maximum  order  specified  was  found  to  be  too  large. 

-3  Corrector  convergence  could  not  be  achieved  for  H . GT.  HMIN. 

-4  Ihe  requested  error  is  smaller  than  can  be  handled  for  this  problem. 

JSTART  An  input  indicator  with  the  following  meanings: 

-1  Repeat  the  last  step  with  a new  H. 

0 Perform  the  first  step.  The  first  step  must  be  done  with  this  value  of 
JSTART  so  that  the  subroutine  can  initialize  itself. 

+ 1 Take  a new  step  continuing  from  the  last. 

JSTART  is  set  to  NQ.  the  current  order  of  the  method  at  exit.  NQ  is  also  the  order 
of  the  maximum  derivative  available. 

MAXPFR  The  maximum  derivative  that  should  be  used  in  the  method.  Since  the  order  is  equal 

to  the  highest  derivative  used,  this  restricts  the  order.  It  must  be  less  than  8 or  7 
for  Adams  or  stiff  methods,  »*espectively . 

PW  A blo«  k of  at  least  N*#2  floating  point  locations. 


r 
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where  n is  the  concentration  (mol/cm  ) of  the  absorbing  species.  Values 
21 

for  Zj,  Z^,  and  a flash-lamp  profile  I(t)  must  be  provided.  is  the 

peak  flash-lamp  output  in  moles  per  unit  volume  of  gas.  Z ^ is  defined  by 
the  relation 


Z 


2 “ 


a l 

v 


n 

a 


3 

cm 

mole 


(14) 


where  a ^ is  the  absorption  coefficient  of  species  n&  in  a gas  container  of 

mean  geometric  length  l.  uj  (WPHOTO)  is  the  average  energy  of  the  ab- 

-1  v 

sorbed  photon  in  cm 

Z j,  Z^.  and  WPHOTO  are  input  by  means  of  the  data  stream  FLASH  1, 
and  the  flash-lamp  output  intensity  profile  is  input  by  the  data  stream  FLASH2. 
The  format  for  the  flash-lamp  intensity  profile  is  as  follows: 


$FLASH2  FLASH  = 1(1),  t(l),  1(2),  t(2) 1(20),  t(20) 


where  I(k)  is  the  flash  output  intensity  at  time  t(k).  I(k)  is  nondimensional, 
and  its  peak  value  is  normalized  to  unity. 

As  shown  in  Appendix  B,  for  an  optically  thin  gas  mixture,  the  initiation 
level  through  photodissociation  of  F^  reduces  to  the  approximate  expression 


— 77  7 

fiyf-  ^ziz2 


(15) 


For  such  cases,  can  be  set  equal  to  one-half  to  obtain 


(16) 
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This  approximation  provides  an  alternative  technique  for  simulating  flash- 
lamp  initiation,  which,  under  certain  conditions,  may  be  easier  to  apply. 

As  an  example,  suppose  parametric  calculations  must  be  made  in  which  it 
is  desired  to  vary  the  flash-lamp  profile  but  maintain  the  same  initiation 
level.  With  this  technique,  after  setting  Z ^ to  the  desired  initiation  level, 
and  Z ^ to  one -half,  all  that  is  needed  for  each  new  run  is  a new  flash -lamp 
profile.  This  procedure  is  certainly  easier  than  using  Eq.  (13)  directly. 

It  should  be  remembered,  however,  that  in  applying  the  above  approximation 
the  flash  profile  is  normalized  in  a different  sense.  In  the  "exact"  treatment 
i.e.,  Eq.  (13),  the  peak  flash  value  is  set  to  unity,  whereas,  in  the  opti- 
cally thin  approximation  case,  the  profile  is  normalized  by 


/l(t)  dt  = 1 

0 


(17) 


V.  CHANGING  RATE  EQUATIONS  AND  COEFFICIENTS 


Substantial  savings  in  computation  time  and  computer  core  storage 
requirements  were  obtained  with  the  present  code  formulation,  wherein  the 
rate  equations  were  input  as  an  integral  part  of  the  program  and  the  rate  co- 
efficients were  precomputed  and  stored  on  disk  for  retrieval  during  program 
execution.  The  drawback  of  this  economy  measure  is  that  changes  in  rate 
equations  or  rate  coefficients  can  require  a rather  significant  effort.  The 
economy  gained,  however,  is  believed  to  justify  this  inconvenience.  The  pro- 
cedure for  making  these  changes  is  described  in  this  section. 

The  basic  steps  involved  in  making  these  changes  follow: 

1.  New  variables  must  be  defined  and  initialized.  The 
total  number  of  integration  variables  is  represented 
by  the  variable  NVAR,  which  should  be  adjusted  to 
reflect  this  change. 

Z.  New  reactions  must  be  incorporated  into  subroutine  DIFFUN, 
and  the  derivatives  affected  by  these  reactions  updated  to 
reflect  the  changes . 

3.  The  thermodynamic  properties  of  any  new  species  intro- 
duced must  be  incorporated  in  the  input  data.  The  new 
species  must  also  be  considered  in  the  gas  temperature 
calculations. 

4.  Pressure -broadening  effects  contributed  by  the  new  species 
may  not  be  negligible,  in  which  case,  the  calculation  for 
the  Lorentz  half -half  width  must  be  revised. 

5.  Rate  modifications  are  made  in  the  rate  data  files.  Rate 
coefficients  corresponding  to  new  reactions  must  be  added 
to  these  data  files. 

6.  Format  can  be  modified  if  desired. 

A detailed  description  is  given  for  the  addition  of  oxygen -related  kinetics 
to  the  code  to  illustrate  the  procedure.  The  reactions  and  rate  constants  uti- 
lized are  essentially  those  compiled  by  Taylor  et  al.  (Ref.  9)  and  given  in 
Table  6.  The  deactivation  of  HF(v)  by  0?  is  taken  to  be  similar  to  that  by  N7, 

L. 

reduced  by  a factor  of  3 (Ref.  10). 
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A. 


NEW  VARIABLES 


Several  new  species  are  introduced  in  these  reactions:  O,  O^. 

OH,  HO^,  and  As  indicated  in  Section  III,  the  variables  in  this  program 

are  represented  by  elements  in  any  array  Y(1,N),  N = 1,  ....  100.  At  this 
point,  the  locations  corresponding  to  N = 15,  16,  and  80  through  100  are  un- 
reserved. Therefore,  the  new  variables  may  be  assigned  to  any  of  these 
positions.  In  this  example,  the  concentrations  of  these  species  will  be  repre- 
sented as  follows: 


Y(l,  15)  = [O] 
Y(l,  16)  = [Oz] 
Y(  1,  80)  = [OH] 
Y(l,  81 ) = [hzO] 
Y(i,  82)  = [F02] 
Y(  1,  83)  =[H02] 


(18) 


For  more  efficient  execution,  it  is  best  to  keep  the  array  dimension  to 
a minimum.  Thus,  whereas  a value  of  N up  to  100  is  possible,  the  present 
example  was  selected  to  occupy  only  up  to  N = 83.  A smaller  array  may  then 
be  utilized  for  the  program  calculations.  The  input  variable  NVAR  was  im- 
plemented to  help  accomplish  this  purpose.  NVAR  designates  the  maximum 
size  of  the  variable  array  that  is  to  be  used  in  the  accompanying  computation. 
Therefore,  in  the  present  situation. 


NVAR  = 83  (19) 

The  new  variables  must  also  be  initialized,  i.e.,  assigned  an  initial  value. 

A zero  value  will  be  automatically  assigned  unless  otherwise  instructed.  It 
is  assumed  for  the  present  purposes  that  only  Y(l,  16)  = [c>2]  has  a nonzero 
value. 
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The  following  procedure  may  be  used  to  introduce  as  a new  input 
variable. 

1.  Expand  the  NAMELIST  group  SPECIES  in  the  subroutine 
INPUT1  to  include  the  variable  02. 

2.  Add  02  to  GASSUM  in  the  subroutine  INPUT  1.  The  frac- 
tional value  of  02  is  then  computed  by  the  statement: 

02  = 02/GASSUM 

3.  Add  02  to  the  common  block  GAS1,  found  in  both  the  main 
program  and  subroutine  INPUT1.  Note  that  in  this  example, 
the  variable  name  02  was  changed  to  R02  in  the  main  program. 

3 

4.  Calculate  the  initial  CL  concentration  (in  moles /cm  ) in  the 
main  program  by  the  statement 

Y(i,  16)  = R02*FCTR  (20) 

where  FCTR  represents  the  appropriate  conversion  factor. 

The  location  of  these  new  statements,  and  others  to  be  described  below, 
are  given  in  the  listings  in  Appendix  C.  The  changes  described  here  are  indi- 
cated in  the  listings  by  the  symbol  O^- 

B.  DIFFUN  MODIFICATIONS 

The  no-oxygen  version  of  the  code  incorporates  150  kinetic  reactions, 
designated  by  the  reaction  numbers  shown  in  Table  A-l.  The  additional  14  re- 
actions considered  here  will  be  designated  by  numbers  151  through  164,  as 
indicated  in  Table  6.  The  variable  NEQN  (main  program)  should  be  set  to  164 
to  reflect  the  increased  number  of  reactions.  For  each  of  these  reactions 
(K),  the  forward  [RCTF(K)|  and  backward  rates  (RCTB(K)j  must  be  com- 
puted. As  an  example,  Reaction  No.  151  is  written  as 

RC  TF  (151)  = ( Y ( 1 , 15)**2)*M12*KFR(151) 

RC  T B ( 1 5 1 ) = M12*Y(1,  16)*KBR(1  51 ) (21) 

where  it  is  recalled  that  Y(l,  15)  and  Y(l,  16)  refer  to  the  species  concentrations 
of  O and  O^,  -espectively . KFR(151)  and  KBR(151)  correspond  to  the  forward 
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and  backward  rate  coefficients,  respectively,  of  Reaction  No.  151,  and  M12 
represents  the  "concentration"  of  the  catalytic  species  O and  O^.  In  the  defi- 
nition of  M12,  the  concentration  of  O was  multiplied  by  2.  82  to  account  for 
the  larger  reaction  cross  section  with  O as  the  third  body  compared  to  02» 
Note  that  M12,  along  with  other  catalytic  species,  is  defined  at  the  beginning 
of  this  subroutine,  and  that  all  these  "M"  variables  must  be  converted  from 
fixed  to  floating  variables  by  the  "REAL"  statement. 

The  subroutine  DIFFUN  will  compute  the  "net  forward  rate"  RF(151) 

RF(151)  = RCTF(151)  - RCTB(151)  (22) 

which  is  then  used  to  calculate  the  contributions  to  the  derivatives  of  each 
species  involved  in  these  reactions.  The  contributions  of  Reaction  Nos.  1 
through  150  have  already  been  accounted  for  in  the  no-oxygen  formulation. 

The  contributions  of  Reaction  Nos.  151  through  164  remain  to  be  included. 

The  reader  can  easily  verify  that  these  reactions  will  affect  the  concentrations 
of  the  following  species: 


[O]  = Y(l,  15) 

[Oz]  = Y (1,  16) 

M = Y(l,  17) 

[H2]  = Y(l,  18) 

[f]  = Y(l,  21)  (23) 

[OH]  = Y(  1,  80) 

[H20]  = Y(l,  81) 

LF02]  = Y(  1,  82) 

[H02]  = Y(l,  83) 


-37- 


Consider  as  an  example  atomic  oxygen  Y(l,  15).  Its  concentration  depends 
directly  on  Reaction  Nos.  151,  152,  154,  155,  158,  159,  160,  and  163.  Thus, 
the  derivative  of  Y(l,  15)  is  written 


DERV1  Y(  15)  = DERV  1Y(  15)  - 2.*[RF(151)  + RF(152)] 

+ RF(  154)  - RF(  1 55 ) - RF(158)  - RF(159)  - RF(160)  (24) 

- RF(  163) 


where  the  first  term  accounts  for  contributions  from  Reaction  Nos.  1 through 
15  0,  normally  computed  in  preceding  sections  of  the  subroutine,  and  subse- 
quent terms  reflect  the  contributions  from  Reaction  Nos.  151  through  164. 
Similar  equations  can  be  written  for  the  other  species.  Two  comments  should 
be  made  in  regard  to  Eq.  (24).  First,  the  leading  term  on  the  right-hand 
side  of  this  equation  is  not  necessary  in  this  instance  since  its  value  is  zero; 
however,  this  situation  is  not  universal,  e.g. , in  the  case  of  DERV1Y(21). 
Thus,  this  term  is  always  suggested  to  reduce  the  possibility  of  errors  result- 
ing from  oversights.  Second,  it  is  possible  to  gain  some  computational  econ- 
omy by  performing  many  of  the  calculations  in  this  equation  separately,  as 
portions  are  repeated  in  subsequent  statements.  But  the  economy  gained  is 
at  the  expense  of  clarity  of  the  individual  rate  equations  (derivatives)  and 
further  complicates  the  process  involved  in  changing  reactions.  The  choice 
made  here  is  to  preserve  that  clarity. 

C.  ENERGY  EQUATION 

The  contributions  of  the  added  species  to  the  terms  of  the  energy  Eq.  (9) 
are  discussed.  The  following  program  variables  are  used. 


CVSUM  = 


N.C 

l 


V. 

l 


where  hh,  Cy.,  and  are,  respectively,  the  concentration,  molar  specific 
heat  at  constant  volume,  and  molar  enthalpy  of  species  i. 

The  values  of  Cy.  and  FL  for  each  of  the  new  species  must  be  incorpo- 
rated into  the  program  data  files  (Section  VI).  The  following  statements  may 
be  used  to  update  the  values  of  the  variables  CVSUM  and  ETHLPY: 

CVSUM  = CVSUM  + CVO(IT)*Y(i,  15)  + CV02(IT)*Y(1,  16) 

+ CVOH(IT)*Y(l,  80)  + CVH20(IT)*Y(1,  81) 

+ CVF02(IT)*Y(1,  82)  + C VH02(IT)*Y(1,  83) 

ETHLPY  = ETHLPY  + EHPYO(IT)*DERVlY(15) 

+ EPHY02(IT)*DERV1Y(16) 

+ EHPYOH(IT)*DERV1Y(80) 

+ EHPYH20(IT)*DERV1Y(81) 

+ EHPYF02(IT)*DERV1Y(82) 

+ EHPYH02(IT)*DERV1  Y(83) 

For  dilute  mixtures  where  concentrations  are  small  (e.g.,  ~0.  5%), 
the  effect  of  these  species  on  the  energy  equation,  as  well  as  on  pressure 
broadening  (as  discussed  in  the  following  paragraphs),  is  not  appreciable. 
Under  these  circumstances,  the  modifications  suggested  in  steps  C and  D may 
be  neglected  with  very  little  error. 

D.  PRESSURE  BROADENING 

The  foreign-gas  broadened  linewidth  is  calculated  in  subroutine  GAIN. 
The  Lorentz  half -half  width  is  represented  by  the  variable  LONTZ.  Hence, 
the  concentration  and  pressure-broadening  coefficients  of  the  new  species 
must  be  incorporated  into  the  computation  of  the  value  of  this  variable.  For 
the  present,  only  will  be  considered. 


The  collision -broadening  characteristics  of  HF  vibrational -rotational 
transitions  by  O.,  are  not  available.  However,  from  published  HC1  pressure- 
broadened  linewidth  data  (Ref.  11),  it  may  be  justifiably  assumed  that  the  O., 
characteristic  is  similar  to  that  of  or  D^.  The  collision  broadening  of 
HF  lines  by  O,  is  taken  to  be  the  same  as  H?  in  the  present  example.  Thus, 
the  following  substitution  is  made  in  the  statement  that  computes  LONTZ: 

H2SUM  — ♦ H2SUM  + Y(l,  16) 

where  H2SUM  is  the  total  H^  concentration  and  Y(i,  16)  is  the  concentration. 
This  modification  will  result  in  a calculated  Lorentz  half-half  width  that  ac- 
counts for  the  contribution. 

E.  RATE-COEFFICIENT  MODIFICATION 

As  stated  previously,  the  rate  coefficients  KFR  and  KBR  used  in  equa- 
tions such  as  Eq.  (21)  were  precomputed  and  stored,  either  on  a magnetic 
tape  or  on  a disk,  or  both.  The  rate  coefficients  are  taken  to  be  of  Arrhenius 
form,  i.  e. , 


Bfr  Efr/RT 
KFR  = A,  T e 
fr 


B,  E,  /RT 

KBR  = A,  T br  e br 
br 


(25) 


KFR  and  KBR  are  connected  by  the  relation 


K - (RT)  r MR 
r (RT)  KBR 


(26) 
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where  is  the  equilibrium  constant  of  the  reaction,  and 


o = V''  v . = (p  . - a . ) 
r Z— i ri  t ri  n 


where  a . and  p . are  stoichiometric  coefficients  of  reaction  r,  as  defined  in 
n ri 

Eq.  (1).  The  equilibrium  constant  is  a function  of  temperature  only  and 
is  given  by 


In  K = - Vv  .H.  + 4-  Y'v  .S. 
r RT  t—i  ri  l R ri  i 


where  FL  (cal/mole)  and  (cal/mole-K)  are,  respectively,  the  molar  enthalpy 
and  molar  entropy  of  species  i.  The  values  of  FF  and  may  be  obtained  from 
the  JANAF  thermodynamic  data  compilation  (Ref.  12). 

The  generation  of  the  rate  coefficient  data  files  is  accomplished  through 
a separate  program,  RATE.  (Note  that  this  is  not  the  subroutine  RATE.  ) 

This  program  computes,  from  input  data,  the  forward  and  backward  rate  co- 
efficients at  25-K  intervals  from  100  to  6000  K.  The  input  data  consist  of  one 
card  for  each  reaction,  which  specifies  the  values  of  A^,  B^,  E^,  A^r»  B^, 
and  E^r  in  the  order  given.  The  format  used  is 

FORMAT  (E15.8,  F9.5,  E16.8,  E15.8,  F9.5,  E16.8) 


The  rate  coefficients  are  calculated  in  this  program  using  Eq.  (25).  The 
generated  data  are  copied  onto  the  logical  files  TAPE7  (KBR)  and  TAPE12  (KFR), 
which  are  subsequently  saved  for  future  use  by  the  model  SPIKE. 

Under  certain  conditions,  a simpler  procedure  may  be  followed  to  ob- 
tain the  needed  rate  coefficients.  For  the  present  example,  the  deactivation 
rate  of  HF(v)  by  was  taken  to  be  similar  to  that  by  N^,  except  that  it  was 
reduced  by  a factor  of  3.  Thus,  the  simulation  of  deactivation  can  be 
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easily  accomplished  by  adding  one -third  of  the  concentration  to  the 
reactions  that  govern  HF(v)  deactivation  by  N.,,  which  is  accomplished  in  the 
subroutine  DIFFUN  by  adjusting  the  definition  of  Mil  to: 

Mil  = SF6  + N2  + 0.3*Y(1,  16) 

where  Y(l,  16)  is  the  element  in  the  variable  array  that  corresponds  to  the 
concentration.  The  rate  coefficients  for  HF  deactivation  by  are  presently 
computed  in  subroutine  RATE.  Thus,  to  ensure  that  this  subroutine  is 
exercised,  the  following  statement  should  be  inserted  in  the  main  program: 

IF(R02.  GT.  0.  ) CALL  RATE  (Y(l,  23)) 
where  R02  refers  to  the  ratio  in  the  initial  gas  mixture. 


VI.  DESCRIPTION  OF  DATA  FILES 


In  addition  to  the  rate  coefficients  contained  in  data  files  TAPE 7 and 
TAPE12,  the  program  must  have  access  to  spectroscopic  and  thermodynamic 
data  that  relate  to  the  species  in  the  reacting  gas  mixture.  These  are  sum- 
marized in  Table  7.  The  data  are  punched  into  data  cards  following  the  for- 
mats indicated  in  the  table  and  input  by  means  of  program  DATA,  which 
reads  the  data  and  stores  them  at  predetermined  locations  (indicated  by  the 
tape  numbers ). 

Most  of  the  data  files  in  the  table  need  no  explanation.  The  rotational 
energies  E(v,  J)  and  transitional  wave  numbers  WC(v,  J)  are  computed  from  the 
data  of  Refs.  7 and  13.  The  rotational  partition  function  can  then  be  computed 
from  the  equation 

Q(v,  T)  = 1 + ]^(2J  + 1)  exp(-  ig-  (29) 

J=i 

The  Einstein  absorption  coefficient  is  computed  from  the  relation 


where  the  matrix  elements  for  the  transitions  ^ ^ are  from  Meredith 

v,  J 

(Ref.  14).  The  heat  capacity  (CV)  and  enthalpy  (EHPY)  data  are  obtained  from 
the  JANAF  thermodynamic  data  compilation. 

The  role  of  pressure  broadening  in  the  pulsed  HF  chemical  laser  is  dis- 
cussed in  detail  in  Ref.  1.  Self -broadened  linewidths  are  calculated  by  the 
use  of  the  data  CCC.  Foreign-gas  broadening  is  based  on  the  experimental 
data  compiled  in  Ref.  1 and  stored  in  array  FB.  The  coefficients  contained 
in  the  array  CCC  were  calculated  by  Meredith  (Ref.  15)  on  the  basis  of  the 
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Anderson  theory  for  pressure  broadening.  The  array  is  dimensioned 
CCC(126,  10),  representing  sets  of  126  coefficients  for  each  of  10  vibrational 
bands.  These  bands  are:  HF(0  - 1),  HF(1  - 2),  HF(2  - 3),  HF(3  - 4), 

HF(4  - 5),  DF(0  - 1),  DF(  1 - 2),  DF(2  - 3),  DF(3  - 4),  and  DF(4  - 5).  Each 
set  of  126  coefficients  is  additionally  partitioned  by  means  of  an  EQUIVALENCE 
statement  into  an  array  of  the  form  C(7,  6,  3).  The  first  index  represents  the 
seven  coefficients  needed  for  the  equation 


/ \ , -m/4  , -m/4  2 -m/2 

•y(m)  = + c^  e + c^me  + c^m  e 

-m^/8  2 -m^/16  -m^/8 

+ Carrie  + c^m  e + c^  e 


(31) 


which  is  the  expression  used  by  Meredith  to  represent  the  results  of  his 
Anderson  theory  calculation.  In  Eq.  (31),  m corresponds  to  the  rotational 
quantum  number  for  the  lower  level  of  the  transition  and  ■y(m)  is  the  Lorentz 
half  linewidth  at  half  maximum  (cm  ^atm  *).  The  coefficients  c.  depend  on 
the  nature  of  the  perturbing  species,  as  well  as  on  the  gas  temperature.  The 
second  index  of  array  C(7,  6,  3)  designates  six  different  possible  perturbers: 
HF(0),  HF ( 1 ),  HF(2),  DF(0),  LF(1),  and  DF(2).  The  third  index  represents 
the  three  temperatures  for  which  the  coefficients  are  given:  300,  600,  and 
900  K.  For  calculations  of  linewidth  at  temperatures  other  than  these,  linear 
interpolations  and  extrapolations  are  presently  used. 

The  foreign-gas  pressure  broadening  coefficients  stored  in  the  array 
FB(15,  8)  correspond  to  the  fifteen  vibrational -rotational  transitions  modeled 
for  each  band  and  eight  possible  perturbing  molecules:  H^,  N^,  F^,  Ar,  He, 

H,  F,  and  SFA.  These  coefficients  were  obtained  experimentally  at  sample 

b 1/2 
gas  temperatures  of  ~300  K.  At  higher  temperatures,  a y(m)  ~ T depen- 
dence is  assumed  for  the  present  calculations,  based  on  a hard  sphere  inter- 
action model  (Ref.  1). 
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APPENDIX  A 


RATE  COEFFICIENTS  FOR  H2  + F2  CHEMICAL  LASER 


The  chemical  kinetic  model  used  for  the  rate -equation  solution  has  been 
suggested  by  Cohen  (Ref.  4)  and  is  shown  in  Table  A-  1.  Rate  coefficients  k and 
k designate  forward  and  backward  rates,  respectively.  For  each  reaction, 
the  missing  rate  coefficient  is  determined  from  the  equilibrium  constant. 


-49- 


Table  A-l.  Rate  Coefficients  for  H,  + F-,  Chemical  Laser 


THIS  PAGE  IS  BEST  QUALITY  PRACTICABLE 
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M)JH  : Vi  + (?)  JH 


Table  A-l.  Rate  Coefficients  for  H + F Chemical  Laser  (Continued) 


e that  dissociation -recombination  reactions  have  been  neglected. 

4.  575  T / 1000  keal/mol 

vibrational  level;  M = collision  partner 
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APPENDIX  B 


FLASH- LAMP  FORMULATION  FOR  OPTICALLY  THIN  GASES 


For  an  optically  thin  gas  mixture,  the  approximation 


exp(-Z2nF2)S1  - Z2nF2 


(B-l) 


can  be  made.  Thus,  Eq.  (29)  can  be  written 


(B-2) 


Since  two  F atoms  are  introduced  for  every  F^  dissociated. 


2 Z^n^Kt) 


(B-3) 


For  low  initiation  levels,  e.g.  , F/F?  S 2%,  Eq.  (B-3)  can  be  integrated 


to  obtain 


mT*ZzlzifW) « 


(B-4) 


where  the  notation  [F^]  was  substituted  for  n^.  . With  a normalized  intensity 
profile  assumed,  i.e.,  ^ 


/l(t)  dt  = 1 


is  obtained.  Thus,  the  quantity  2 Z^Z^  is  approximately  the  initiation  level  of 
the  laser,  and  can  be  empirically  set  to  the  measured  experimental  value. 
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THE  IVAN  A.  GETTING  LABORATORIES 


The  Laboratory  Operations  of  The  Aerospace  Corporation  is  conducting 
experimental  and  theoretical  investigations  necessary  for  the  evaluation  and 
application  of  scientific  advances  to  new  military  concepts  and  systems.  Ver- 
satility and  flexibility  have  been  developed  to  a high  degree  by  the  laboratory 
personnel  in  dealing  with  the  many  problems  encountered  in  the  nation's  rapidly 
developing  space  and  missile  systems.  Expertise  in  the  latest  scientific  devel- 
opments is  vital  to  the  accomplishment  of  tasks  related  to  these  problems.  The 
laboratories  that  contribute  to  this  research  are: 

Aerophysics  Laboratory:  Launch  and  reentry  aerodynamics,  heat  trans- 
fer, reentry  physics,  chemical  kinetics,  structural  mechanics,  flight  dynamics, 
atmospheric  pollution,  and  high-power  gas  lasers. 

Chemistry  and  Physics  Laboratory:  Atmospheric  reactions  and  atmos- 
pheric  optics,  chemical  reactions  in  polluted  atmospheres,  chemical  reactions 
of  excited  species  in  rocket  plumes,  chemical  thermodynamics,  plasma  and 
laser-induced  reactions,  laser  chemistry,  propulsion  chemistry,  space  vacuum 
and  radiation  effects  on  materials,  lubrication  and  surface  phenomena,  photo- 
sensitive materials  and  sensors,  high  precision  laser  ranging,  and  the  appli- 
cation of  physics  and  chemistry  to  problems  of  law  enforcement  and  biomedicine. 

Electronics  Research  Laboratory:  Electromagnetic  theory,  devices,  and 
propagation  phenomena,  including  plasma  electromagnetics;  quantum  electronics, 
lasers,  and  electro-optics;  communication  sciences,  applied  electronics,  semi- 
conducting, superconducting,  and  crystal  device  physics,  optical  and  acoustical 
imaging;  atmospheric  pollution;  millimeter  wave  and  far-infrared  technology. 

Materials  Sciences  Laboratory:  Development  of  new  materials;  metal 
matrix  composites  and  new  forms  of  carbon;  test  and  evaluation  of  graphite 
and  ceramics  in  reentry;  spacecraft  materials  and  electronic  components  in 
nuclear  weapons  environment;  application  of  fracture  mechanics  to  stress  cor- 
rosion and  fatigue-induced  fractures  in  structural  metals. 

Space  Sciences  Laboratory:  Atmospheric  and  ionospheric  physics,  radia- 
tion from  the  atmosphere,  density  and  composition  of  the  atmosphere,  aurorae 
and  airglow;  magnetospheric  physics,  cosmic  rays,  generation  and  propagation 
of  plasma  waves  in  the  magnetosphere;  solar  physics,  studies  of  solar  magnetic 
fields;  space  astronomy,  x-ray  astronomy;  the  effects  of  nuclear  explosions, 
magnetic  storms,  and  solar  activity  on  the  earth's  atmosphere,  ionosphere,  and 
magnetosphere;  the  effects  of  optical,  electromagnetic,  and  particulate  radia- 
tions in  space  on  space  systems. 

THE  AEROSPACE  CORPORATION 
El  Segundo,  California 


» 


